Problem Set 3 Key: General(ized) Linear Models

PLS 206 — Applied Multivariate Modeling

Tip

Try the problems yourself before consulting this key. Every number below was checked in R against data/pistachio.csv (300 trees, 7 columns) on 2026-10-12.

library(ggplot2)
library(emmeans)
pistachio <- read.csv("https://greymonroe.github.io/PLS_206/data/pistachio.csv")
str(pistachio)   # 300 obs. of 7 variables

---

Part 1: Genotype and Yield

Q1. (Intercept) = 10.49, the mean yield of chaparillo (the reference genotype, first in the alphabet). genotypekerman = +0.19: kerman trees yield 0.19 more than chaparillo trees on average (kerman mean 10.68).

ggplot(pistachio, aes(genotype, yield)) + geom_boxplot() + geom_jitter(width = 0.1, alpha = 0.4) + theme_classic()
m1 <- lm(yield ~ genotype, data = pistachio)
summary(m1)
tapply(pistachio$yield, pistachio$genotype, mean)   # 10.49, 10.68, 9.39

Q2. No evidence that the genotypes differ in yield: F = 1.19 on 2 and 297 df, p = 0.31. Tukey: no pair differs (all adjusted p > 0.3).

anova(m1)
TukeyHSD(aov(m1))

---

Part 2: Disease, Yield, and Genotype

Q3. Yield lost to disease: chaparillo 12.76, t41 7.62, kerman 4.50 (all p < 0.0001). Kerman loses the least, 4.5 (13.83 healthy vs 9.33 diseased).

pistachio$disease_f <- factor(pistachio$disease)
ggplot(pistachio, aes(genotype, yield, fill = disease_f)) + geom_boxplot() + theme_classic()
m2 <- lm(yield ~ genotype * disease_f, data = pistachio)
summary(m2)
emmeans(m2, ~ disease_f | genotype)
pairs(emmeans(m2, ~ disease_f | genotype))

Q4. With the interaction in the model, genotypekerman compares kerman with chaparillo among healthy trees only (disease_f = 0 is the reference level): healthy kerman yields 3.80 less than healthy chaparillo (13.83 vs 17.63), a clear difference. The genotype row of anova() tests the genotype term as a whole, entered first in a sequential table: the genotypes compared ignoring disease, which is the Q2 comparison (genotype means of 10.49, 10.68, 9.39) and not significant. The genotype differences reverse between healthy and diseased trees and cancel out on average.

anova(m2)   # genotype row: F = 2.25, p = 0.107

Q5. Yes: the genotype:disease_f interaction is strong (F = 18.5 on 2 and 294 df, p = 2.6 × 10⁻⁸). With the disease present, kerman: diseased kerman trees yield 9.33 vs 4.88 (chaparillo) and 4.51 (t41), so it tolerates the disease best. If the disease can be prevented, chaparillo: healthy chaparillo trees yield the most (17.63 vs 13.83 kerman and 12.13 t41). Any answer that uses the interaction and the cell means this way is correct; noting that kerman is also the genotype most often diseased (Q6) is a plus.

---

Part 3: Disease Rates

Q6. Proportion diseased: chaparillo 0.56, kerman 0.70, t41 0.36. The model reproduces them exactly: plogis(0.2412) = 0.56 (chaparillo), plogis(0.2412 - 0.8165) = 0.36 (t41).

prop.table(table(pistachio$genotype, pistachio$disease), 1)
g <- glm(disease ~ genotype, family = binomial, data = pistachio)
summary(g)
plogis(coef(g)[1])                  # chaparillo: 0.56
plogis(coef(g)[1] + coef(g)[3])     # t41: 0.36

Q7. Odds ratio 0.44, 95% CI 0.25 to 0.78. Example: “T41 trees had lower odds of disease than chaparillo trees (odds ratio 0.44, 95% CI 0.25 to 0.78; z = −2.82, p = 0.005; logistic regression, n = 300).” A sentence giving the odds ratio, its CI, and the test is enough. (T41 has less than half the odds of disease.)

exp(coef(g))     # genotypet41 0.44
exp(confint(g))  # 0.25 to 0.78

---

Part 4: An Early Warning from Drones

Q8. March NDVI: slope −0.81 (p = 0.034), accuracy 0.563. Guessing “diseased” for every tree is right 0.54 of the time. March NDVI is at most weakly associated with disease (p = 0.03) and adds almost nothing as a predictor (56% vs 54%): not useful.

m_mar <- glm(disease ~ NDVI_march, family = binomial, data = pistachio)
summary(m_mar)
pred <- as.numeric(fitted(m_mar) > 0.5)
mean(pred == pistachio$disease)   # 0.563
mean(pistachio$disease == 1)      # 0.54

Q9. Accuracy: March + June 0.743; March + June + August 0.813. In the March + June model the March slope is +9.4 and the June slope is −10.0, nearly equal and opposite: the model flags trees whose greenness dropped from March to June (green in spring, browning by summer). A model on the change alone (NDVI_june − NDVI_march, slope −10.1) is 0.74 accurate, the same as the two-variable model.

m_jun <- glm(disease ~ NDVI_march + NDVI_june, family = binomial, data = pistachio)
summary(m_jun)
mean(as.numeric(fitted(m_jun) > 0.5) == pistachio$disease)   # 0.743

m_aug <- glm(disease ~ NDVI_march + NDVI_june + NDVI_august, family = binomial, data = pistachio)
summary(m_aug)
mean(as.numeric(fitted(m_aug) > 0.5) == pistachio$disease)   # 0.813

pistachio$change <- pistachio$NDVI_june - pistachio$NDVI_march
m_chg <- glm(disease ~ change, family = binomial, data = pistachio)
mean(as.numeric(fitted(m_chg) > 0.5) == pistachio$disease)   # 0.74

Q10. Any reasoned answer. Example: after June. March is no better than guessing; June reaches 74% (20 points above the baseline) with the rest of the season left to treat; August adds 7 points but leaves little time for the treatment to recover yield. If treatment is very expensive and late treatment still helps, waiting for August is defensible.